
## Playing Politics with Environmental Protection: The Political Economy of Designating Protected Areas
## Jorge Mangonnet, Jacob Kopas, and Johannes Urpelainen
## August 2021

## File: code for estimating and storing regression models
## NOTE: run this code before table_figures.R


### Preamble -------------------------------------------------------------------

## Clean workspace
rm(list=ls())		
dev.off()       
cat("\014")     

## Load packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  maptools,
  spatstat,
  sp,
  spdep,
  splm, 
  rgdal, 
  foreign,
  stargazer,
  xtable,
  ggplot2,
  plyr,
  dplyr, 
  reshape2,
  lfe,
  plm,
  stringr, 
  magrittr, 
  rdrobust, 
  rddensity
)

## Set user
jk <- "~/Dropbox/Transamazônica/"
jm <- "~/Dropbox/Research/Transamazonica/"
user <- jm

## Functions for regression models and plots
source(paste0(user, "manuscript/replication/Brazil_func.R"))



### LOAD, RECODE, AND SUBSET DATA -----------------------------------------

## Load main data frame
df <- read.csv(paste0(user, "manuscript/replication/data_main.csv"))

## Geographic
# Distance measure: convert meters to kilometers
df$dist <- df$nb_d / 1000
# Coordinate polynomials
df$xy <- df$x * df$y
df$x2 <- df$x^2
df$y2 <- df$y^2
df$x3 <- df$x^3
df$y3 <- df$y^3
df$x2y <- df$x^2 * df$y
df$xy2 <- df$x * df$y^2
# Bandwidth (for rdrobust function)
df$bw.dist <- ifelse(df$p_coal == 0, df$dist*-1, df$dist)

## Subset: strict model
# Coalition
df.s <- df[which(df$coal_match == 1), ]
# Party
df.p <- df[which(df$party_match == 1), ]

## Comparing 20km, 15km, and 10km bandwidth
df.20 <- df.s[which(df.s$nb_d <= 20000), ]  # 20 km
df.15 <- df.s[which(df.s$nb_d <= 15000), ]  # 15 km
df.10 <- df.s[which(df.s$nb_d <= 10000), ]  # 10 km

## Subset: governor-mayor model (only grid cells w/in states)
# Coalition
df.state.c <- df[which(df$state_coal_match == 1), ]
df.state.c <- subset(df.state.c, substr(df.state.c$codigo, 0, 2) == substr(df.state.c$nb, 0, 2))
# Party
df.state.p <- df[which(df$state_party_match == 1), ] 
df.state.p <- subset(df.state.p, substr(df.state.p$codigo, 0, 2) == substr(df.state.p$nb, 0, 2))

## Subset: president-governor model (only grid cells b/ states)
# Coalition
df.prgo.c <- df[which(df$prgo_coal_match == 1), ]
df.prgo.c <- subset(df.prgo.c, substr(df.prgo.c$codigo, 0, 2) != substr(df.prgo.c$nb, 0, 2))
df.prgo.c$UF <- substr(df.prgo.c$codigo, start = 1, stop = 2)
df.prgo.c$nb_UF <- substr(df.prgo.c$nb, start = 1, stop = 2)
df.prgo.c$s_pair <- paste0(df.prgo.c$UF,"_",df.prgo.c$nb_UF)
# Party
df.prgo.p <- df[which(df$prgo_party_match == 1), ]
df.prgo.p <- subset(df.prgo.p, substr(df.prgo.p$codigo, 0, 2) != substr(df.prgo.p$nb, 0, 2))
df.prgo.p$UF <- substr(df.prgo.p$codigo, start = 1, stop = 2)
df.prgo.p$nb_UF <- substr(df.prgo.p$nb, start = 1, stop = 2)
df.prgo.p$s_pair <- paste0(df.prgo.p$UF,"_",df.prgo.p$nb_UF)

## Subset: presidential party period
# Cardoso (PSDB)
df.card <- df.s[df.s$year %in% 1997:2002, ]
# Lula/Dilma (PT)
df.lula <- df.s[df.s$year %in% 2003:2012, ]

## Subset: drop saturated cell data set (Y = 1.0)
df.sat <- df.s
d <- rep(NA, nrow(df.sat))
pb <- txtProgressBar(1, nrow(df.sat), style = 3)
for (i in 1:nrow(df.sat)){
  setTxtProgressBar(pb, i)
  if (df.sat$federal[i] < 1){
    d[i] <- 0
  } else {
    sel <- which(df.sat$ID == df.sat$ID[i] & df.sat$year > df.sat$year[i])
    d[sel] <- 1
    if (is.na(d[i])){
      d[i] <- 0
    }
  } 
}
df.sat <- df.sat[d == 0, ]

## Subset: by mayoral electoral year (for covariate balance tests)
df_1996 <- subset(df.s, year == 1997)
df_2000 <- subset(df.s, year == 2001)
df_2004 <- subset(df.s, year == 2005)
df_2008 <- subset(df.s, year == 2009)

##  Municipal-level data sets
# Elections
df.muni <- read.csv(paste0(user, "data/analysis/grid_yr dataset/data_muniborder.csv"))
df.muni <- df.muni[which(df.muni$coal_match == 1),]
# Economic extraction
df.ex <- read.csv(paste0(user, "/data/analysis/grid_yr dataset/data_muniextractive.csv"))
df.ex <- subset(df.ex, codigo %in% df.muni$codigo)



### (A) MAIN MODELS ------------------------------------------------------------
### "Findings" section of the paper

## Table 1: Federal Protected Areas and Political Alignment ----
## NOTE: only matching years for treatment varying cells
# Coalition Alignment (Models 1, 2, 3)
formulas.c <- list(federal ~ p_coal + federal_pre97 | m_pair | 0 | m_pair,
                   federal ~ p_coal | m_pair + ID | 0 | m_pair,
                   federal ~ p_coal | m_pair + ID + state_yr | 0 | m_pair)
# Party Alignment (Models 4, 5, 6)
formulas.p <- list(federal ~ p_align + federal_pre97 | m_pair | 0 | m_pair,
                   federal ~ p_align | m_pair + ID | 0 | m_pair,
                   federal ~ p_align | m_pair + ID + state_yr | 0 | m_pair)

## Estimate and store (only matching years for treatment varying cells)
# Coalition Alignment (Models 1, 2, 3)
main.model(formulas = formulas.c, code = "coal", data = df.s)
coal.ug <- rep(length(unique(df.s$ID)), length(coal))
# Party Alignment (Models 4, 5, 6)
main.model(formulas = formulas.p, code = "party", data = df.p)
party.ug <- rep(length(unique(df.p$ID)), length(party))


## Table 2: Presidential Vote Share and Federal Protected Areas ----
formulas.pres <- list(pres_vote ~ federal + federal_pre97 | 0 | 0 | codigo,
                      pres_vote ~ federal | codigo | 0 | codigo,
                      pres_vote ~ federal | codigo + state_yr | 0 | codigo,
                      pres_vote ~ federal + federal_pre97 + p_coal + p_coal:federal | 0 | 0 | codigo,
                      pres_vote ~ federal + p_coal + p_coal:federal | codigo | 0 | codigo,
                      pres_vote ~ federal + p_coal + p_coal:federal | codigo + state_yr | 0 | codigo)
# Estimate and store
main.model(formulas.pres, "pres", df.muni)
pres.ug <- rep(nrow(unique(df.muni[ , c("codigo", "m_pair")])), length(pres))



### (B) EXPLORATORY ANALYSES ---------------------------------------------------
### Appendix Section A6

## Table A5: Comparison of Extensive and Intensive Margins of Coalition Alignment ----
formulas.muni <- list(federal ~ p_coal + federal_pre97 | m_pair | 0 | m_pair,
                      federal ~ p_coal | m_pair + codigo | 0 | m_pair,
                      federal ~ p_coal | m_pair + codigo + state_yr | 0 | m_pair,
                      federal_prop ~ p_coal + federal_prop_pre97 | m_pair | 0 | m_pair,
                      federal_prop ~ p_coal | m_pair + codigo | 0 | m_pair,
                      federal_prop ~ p_coal | m_pair + codigo + state_yr | 0 | m_pair)
# Estimate and store
main.model(formulas.muni, "muni", df.muni)
muni.ug <- rep(nrow(unique(df.muni[ , c("codigo", "m_pair")])), length(muni))


## Table A6: Interaction between Coalition Alignment and Presidential Vote Share ----
## NOTE: Figure 3 in paper based on these models
formulas.cs <- list(federal ~ p_coal + federal_pre97 + p_coal:pres_vote_pre + pres_vote_pre | m_pair | 0 | m_pair,
                    federal ~ p_coal + p_coal:pres_vote_pre + pres_vote_pre | m_pair + ID | 0 | m_pair,
                    federal ~ p_coal + p_coal:pres_vote_pre + pres_vote_pre | m_pair + ID + state_yr | 0 | m_pair)
# Estimate and store
main.model(formulas.cs, "core", df.s)
core.ug <- rep(length(unique(df.s$ID)), length(core))


## 10-point bins of vote share ----
## NOTE: model for Appendix Figure A6 (regression table NOT IN PAPER)
## WARNING: computationally demanding, takes a long time to complete (~1 hour)
formulas.cbin <- list(federal ~ p_coal + p_coal:pres_30 + pres_30 + 
                        p_coal:pres_40 + pres_40 + p_coal:pres_50 + pres_50 +
                        p_coal:pres_60 + pres_60 + p_coal:pres_70 + pres_70 + 
                        p_coal:pres_80 + pres_80 + p_coal:pres_90 + pres_90 + 
                        p_coal:pres_100 + pres_100 | m_pair + ID + state_yr | 0 | m_pair)
# Estimate and store
main.model(formulas.cbin, "cbin", df.s)
cbin.ug <- rep(length(unique(df.s$ID)), length(cbin))


## Table A7: Interaction between Coalition Alignment and Margin of Victory ----
formulas.vm <- list(federal ~ p_coal + federal_pre97 + pres_marg_pre + p_coal:pres_marg_pre  | m_pair | 0 | m_pair,
                    federal ~ p_coal + pres_marg_pre + p_coal:pres_marg_pre  | m_pair + ID | 0 | m_pair,
                    federal ~ p_coal + pres_marg_pre + p_coal:pres_marg_pre  | m_pair + ID + state_yr | 0 | m_pair)
# Estimate and store
main.model(formulas.vm, "vmarg", df.s)
vmarg.ug <- rep(length(unique(df.s$ID)), length(vmarg))


## Table A8: Difference-in-Differences Estimation of the Effect of Protected Areas on Local Extractive Industries ----
## NOTE: Figure 4 based on these models
formulas.ex <- list(log(soybeantn + 1) ~ protected_prop + post2001 + I(protected_prop * post2001) | 0 | 0 | codigo,
                    log(soybeantn + 1) ~ post2001 + I(protected_prop * post2001) | codigo | 0 | codigo,
                    log(soybeantn + 1) ~ I(protected_prop * post2001) | codigo + state_yr | 0 | codigo,
                    log(minleases_den + 1) ~ protected_prop + post2001 + I(protected_prop * post2001) | 0 | 0 | codigo,
                    log(minleases_den + 1) ~ post2001 + I(protected_prop * post2001) | codigo | 0 | codigo,
                    log(minleases_den + 1) ~ I(protected_prop * post2001) | codigo + state_yr | 0 | codigo)
# Estimate and store
main.model(formulas.ex, "ex", df.ex)
ex.ug <- rep(nrow(unique(df.ex[ , c("codigo")])), length(ex))


## Table A9: Interaction between Coalition Alignment and Potential for Economic Exploitation ----
formulas.int <- list(federal ~ p_coal + federal_pre97 +
                       p_coal:deforested + deforested | m_pair | 0 | m_pair,
                     federal ~ p_coal + federal_pre97 +
                       p_coal:pastures + pastures | m_pair | 0 | m_pair,
                     federal ~ p_coal + federal_pre97 + 
                       p_coal:soybeans + soybeans | m_pair | 0 | m_pair)
# Estimate and store
main.model(formulas.int, "int", df.s)
int.ug <- rep(length(unique(df.s$ID)), length(int))


## Table A10: Interaction between Coalition Alignment and Federal Roadways ----
formulas.int.rds <- list(federal ~ p_coal + federal_pre97 + 
                           p_coal:br230_distance_1993 + br230_distance_1993 | m_pair | 0 | m_pair,
                         federal ~ p_coal + federal_pre97 + 
                           p_coal:main_rds_dist + main_rds_dist | m_pair | 0 | m_pair, 
                         federal ~ p_coal + federal_pre97 + 
                           p_coal:rds_dist + rds_dist | m_pair | 0 | m_pair)
# Estimate and store
main.model(formulas.int.rds, "int.rds", df.s)
int.rds.ug <- rep(length(unique(df.s$ID)), length(int.rds))



### (C) PLACEBO TESTS AND ROBUSTNESS CHECKS ------------------------------------
### Appendix Section A7

## Table A11: Effect of Coalition Alignment on Federal Protected Areas, State Protected Areas, and Indigenous Lands ----
## NOTE: Figure 5 based on these models
## WARNING: computationally demanding, takes a long time to complete (~1 hour)
formulas.type <- list(federal ~ p_coal + federal_pre97 | m_pair | 0 | m_pair,
                      federal ~ p_coal | m_pair + ID | 0 | m_pair,
                      federal ~ p_coal | m_pair + ID + state_yr | 0 | m_pair,
                      indigen ~ p_coal + indigen_pre97 | m_pair | 0 | m_pair,
                      indigen ~ p_coal | m_pair + ID | 0 | m_pair,
                      indigen ~ p_coal | m_pair + ID + state_yr | 0 | m_pair,
                      state ~ p_coal + state_pre97 | m_pair | 0 | m_pair,
                      state ~ p_coal | m_pair + ID | 0 | m_pair,
                      state ~ p_coal | m_pair + ID + state_yr | 0 | m_pair)
# Estimate and store
main.model(formulas.type, "type", df.s)
type.ug <- rep(length(unique(df.s$ID)), length(type))


## Table A12: State Protected Areas and Governor-Mayor Alignment ----
# Coalition Alignment (Models 1, 2, 3)
formulas.sc <- list(state ~ s_coal + state_pre97 | m_pair | 0 | m_pair,
                    state ~ s_coal | m_pair + ID | 0 | m_pair,
                    state ~ s_coal | m_pair + ID + state_yr | 0 | m_pair)
# Party Alignment (Models 4, 5, 6)
formulas.sp <- list(state ~ s_align + state_pre97 | m_pair | 0 | m_pair,
                    state ~ s_align | m_pair + ID | 0 | m_pair,
                    state ~ s_align | m_pair + ID + state_yr | 0 | m_pair)

## Estimate and store
# Coalition Alignment (Models 1, 2, 3)
main.model(formulas.sc, "s_coal", df.state.c)
scoal.ug <- rep(length(unique(df.state.c$ID)), length(s_coal))
# Party Alignment (Models 4, 5, 6)
main.model(formulas.sp, "s_party", df.state.p)
sparty.ug <- rep(length(unique(df.state.p$ID)), length(s_party))


## Table A13: Controlling for Distance to the Border ----
formulas.rd <- list(federal ~ p_coal + federal_pre97 + dist | m_pair + codigo + state_yr | 0 | m_pair, 
                    federal ~ p_coal + federal_pre97 + dist + dist:p_coal | m_pair + codigo + state_yr | 0 | m_pair)
# Estimate and store
main.model(formulas.rd, "rd", df.s)
cxrd.ug <- rep(length(unique(df.s$ID)), length(rd))


## Table A14: Controlling for Geographic Coordinates ----
## NOTE: Dell's (2010) semiparametric approach
formulas.dell <- list(federal ~ p_coal + federal_pre97 | m_pair + state_yr | 0 | m_pair,
                      federal ~ p_coal + federal_pre97 + x + y | m_pair + state_yr | 0 | m_pair,
                      federal ~ p_coal + federal_pre97 + x + y + x2 + y2 + xy | m_pair + state_yr | 0 | m_pair,
                      federal ~ p_coal + federal_pre97 + x + y + x2 + y2 + xy + x3 + y3 + x2y + xy2 | m_pair + state_yr | 0 | m_pair)
## Estimate and store
main.model(formulas.dell, "dell", df.s)
dell.ug <- rep(length(unique(df.s$ID)), length(dell))


## Table A15: Nonparametric Estimation of the Effect of Coalition Alignment ----
## NOTE (a): Keele & Titiunik's (2015) nonparametric approach
## NOTE (b): estimates are directly stored in the tables_figures.R script
## WARNING: computationally demanding
np <- rdrobust(y = df.s$federal,
               x = df.s$bw.dist, 
               h = 25, 
               c = 0, 
               masspoints = 'adjust')
# Check coeffs, SEs, and CIs 
np[["coef"]]
np[["se"]]
np[["ci"]]


## Table A16: Comparison of Main Models at Different Bandwidths ----
## NOTE: Full model (all three models) / same formula as main models
# Estimate and store
sub.model(formulas.c, "bands", list(df.20, df.15, df.10))
bands.ug <- c(rep(length(unique(df.20$ID)), 3),
              rep(length(unique(df.15$ID)), 3),
              rep(length(unique(df.10$ID)), 3))


## Table A17: Main Models with Full Panel of Grid Cell-Year Observations ----
## NOTE: all years with at least one matching year
formulas.cfull <- list(federal ~ p_coal + federal_pre97 | m_pair + year | 0 | m_pair,
                       federal ~ p_coal | m_pair + ID + year | 0 | m_pair,
                       federal ~ p_coal | m_pair + ID + state_yr | 0 | m_pair)
# Estimate and store
main.model(formulas = formulas.cfull, code = "full", data = df)
full.ug <- rep(length(unique(df$ID)), length(full))


## Table A18: Main Models without Fully Saturated Grid Cells ----
## NOTE: same formula as main models
# Estimate and store
main.model(formulas = formulas.c, code = "sat", data = df.sat)
sat.ug <- rep(length(unique(df.sat$ID)), length(sat))


## Table A19: Main Models with Imbalanced Pre-treatment Covariates as Control Variables ----
## NOTE: included covariates that failed balance test at p<.10
formulas.ccv <- list(federal ~ p_coal + federal_pre97 + ideo_imp_94 +
                       deforested + urb_extent + sgi_mean + elf + 
                       cfr_mean + acc_mean + rcr_mean + sbi_mean + 
                       sbr_mean + mml_mean + cfi_mean + rci_mean | m_pair | 0 | m_pair,
                     federal ~ p_coal + federal_pre97 + ideo_imp_94 +
                       cfr_mean + acc_mean + rcr_mean + sbi_mean + 
                       sbr_mean + mml_mean + cfi_mean + rci_mean | m_pair + state_yr | 0 | m_pair)
# Estimate and store
main.model(formulas = formulas.ccv, code = "ccv", data = df.s)
ccv.ug <- rep(length(unique(df.s$ID)), length(ccv))



### (D) EXTENSIONS -------------------------------------------------------------
### Appendix Section A8

## Table A20: Environmental Embargoes and Political Alignment ----
# Coalition Alignment
formulas.embc <- list(embargoed.area ~ p_coal + federal_pre97 | m_pair | 0 | m_pair,
                      embargoed.area ~ p_coal | m_pair + ID | 0 | m_pair,
                      embargoed.area ~ p_coal | m_pair + ID + state_yr | 0 | m_pair)
# Party Alignment
formulas.embp <- list(embargoed.area ~ p_align + federal_pre97 | m_pair | 0 | m_pair,
                      embargoed.area ~ p_align | m_pair + ID | 0 | m_pair,
                      embargoed.area ~ p_align | m_pair + ID + state_yr | 0 | m_pair)

## Estimate and store
# Coalition Alignment
main.model(formulas.embc, "embc", df.s)
embc.ug <- rep(length(unique(df.s$ID)), length(embc))
# Party Alignment
main.model(formulas.embp, "embp", df.p)
embp.ug <- rep(length(unique(df.p$ID)), length(embp))


## Table A21: Federal Protected Areas and Coalition Alignment in the PSDB Presidential Term (1997-2002) ----
## NOTE: same formula as in main models
# Estimate and store
main.model(formulas.c, "cardoso", df.card)
card.ug <- rep(length(unique(df.card$ID)), length(cardoso))


## Table A22: Federal Protected Areas and Coalition Alignment in the PT Presidential Terms (2003-2012) ----
## NOTE: same formula as in main models
# Estimate and store
main.model(formulas.c, "lula", df.lula)
lula.ug <- rep(length(unique(df.lula$ID)), length(lula))


## Table A23: Incumbent Mayor’s Vote Share and Federal Protected Areas ----
formulas.mayor <- list(mayor_vote ~ federal + federal_pre97 | 0 | 0 | codigo,
                       mayor_vote ~ federal | codigo | 0 | codigo,
                       mayor_vote ~ federal | codigo + state_yr | 0 | codigo,
                       mayor_vote ~ federal + federal_pre97 + p_coal + p_coal:federal | 0 | 0 | codigo,
                       mayor_vote ~ federal + p_coal + p_coal:federal | codigo | 0 | codigo,
                       mayor_vote ~ federal + p_coal + p_coal:federal | codigo + state_yr | 0 | codigo)
# Estimate and store  
main.model(formulas.mayor, "mayor", df.muni)
mayor.ug <- rep(nrow(unique(df.muni[ , c("codigo", "m_pair")])), length(mayor))


## Table A24: Difference-in-Differences Estimation of the Effect of Protected Areas on Local Peasant Agriculture ----
formulas.pea <- list(log(familiesn + 1) ~ protected_prop + post2001 + I(protected_prop * post2001) | 0 | 0 | codigo,
                     log(familiesn + 1) ~ post2001 + I(protected_prop * post2001) | codigo | 0 | codigo,
                     log(familiesn + 1) ~ I(protected_prop * post2001) | codigo + state_yr | 0 | codigo, 
                     log(peasantf + 1) ~ protected_prop + post2001 + I(protected_prop * post2001) | 0 | 0 | codigo,
                     log(peasantf + 1) ~ post2001 + I(protected_prop * post2001) | codigo | 0 | codigo,
                     log(peasantf + 1) ~ I(protected_prop * post2001) | codigo + state_yr | 0 | codigo)
# Estimate and store
main.model(formulas.pea, "pea", df.ex)
pea.ug <- rep(nrow(unique(df.ex[ , c("codigo")])), length(pea))


## Table A25: Interaction between Coalition Alignment and Municipal Deforestation ----
formulas.md <- list(federal ~ p_coal + federal_pre97 + deforested_muni_pr + p_coal:deforested_muni_pr | m_pair | 0 | m_pair, 
                    federal ~ p_coal + deforested_muni_pr + p_coal:deforested_muni_pr | m_pair + ID | 0 | m_pair, 
                    federal ~ p_coal + deforested_muni_pr + p_coal:deforested_muni_pr | m_pair + ID + state_yr | 0 | m_pair, 
                    federal ~ p_coal + federal_pre97 + deforested_muni_ca + p_coal:deforested_muni_ca | m_pair | 0 | m_pair, 
                    federal ~ p_coal + deforested_muni_ca + p_coal:deforested_muni_ca | m_pair + ID | 0 | m_pair, 
                    federal ~ p_coal + deforested_muni_ca + p_coal:deforested_muni_ca | m_pair + ID + state_yr | 0 | m_pair)
# Estimate and store
main.model(formulas.md, "md", df.s)
md.ug <- rep(length(unique(df.s$ID)), length(md))


## Table A26: Federal Protected Areas and President-Governor Alignment ----
formulas.gc <- list(federal ~ g_coal + federal_pre97 | s_pair | 0 | s_pair,
                    federal ~ g_coal | s_pair + ID | 0 | s_pair,
                    federal ~ g_coal | s_pair + ID + year | 0 | s_pair)
# Estimate and store
main.model(formulas.gc, "g_coal", df.prgo.c)
gcoal.ug <- rep(length(unique(df.prgo.c$ID)), length(g_coal))



### (E) IDENTIFYING ASSUMPTIONS ------------------------------------------------
### Regressions for balance plots and tables

## Figure A3: Balance test for Coalition Alignment on 40 pre-treatment covariates in the 1996 mayoral election ----
## NOTE: regress treatment on each pre-treatment covariates by mayoral electoral year
c96.ideol <- felm(ideo_imp ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.defo <- felm(deforested ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.rain <- felm(rain_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.eva <- felm(scale(eva_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.temp <-felm(temp_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.climaggr <- felm(scale(climate_aggr_lv) ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.alti <- felm(altg_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.slope <- felm(scale(slp_index) ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.access <- felm(scale(acc_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.workab <- felm(scale(wrk_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.nutri <- felm(scale(nut_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.water <- felm(water_area ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.veget <- felm(vegetation_area ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.cacao_i <- felm(cci_mean ~ p_coal| m_pair | 0 | m_pair, data=df_1996)
c96.cacao_r <- felm(ccr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.coffee_i <- felm(cfi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.coffee_r <- felm(cfr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.grass_i <- felm(pgi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.grass_r <- felm(pgr_mean ~ p_coal| m_pair | 0 | m_pair, data=df_1996)
c96.legum_i <- felm(pli_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.legum_r <- felm(plr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.maize_i <- felm(mzi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.maize_r <- felm(mzr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.rice_i <- felm(rci_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.rice_r <- felm(rcr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.soy_i <- felm(sbi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.soy_r <- felm(sbr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.sugar_i <- felm(sgi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.sugar_r <- felm(sgr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.amphib <- felm(amp_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.birds <- felm(brd_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.mammal <- felm(mml_mean ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.urbdist <- felm(HubDist ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.urbexte <- felm(urb_extent ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.br230 <- felm(br230_distance_1993 ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.mainrds <- felm(main_rds_dist ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.roads <- felm(rds_dist ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.elf <- felm(elf ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.popc <- felm(pop95c_avg ~ p_coal | m_pair | 0 | m_pair, data=df_1996)
c96.popd <- felm(pop95d_avg ~ p_coal | m_pair | 0 | m_pair, data=df_1996)

## Table A3: Balance Test for Coalition Alignment on 40 Pre-Treatment Covariates in the 2000, 2004, and 2008 Mayoral Elections ----
# Balance - 2000
c00.ideol <- felm(ideo_imp ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.defo <- felm(deforested ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.rain <- felm(rain_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.eva <- felm(scale(eva_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.temp <-felm(temp_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.climaggr <- felm(scale(climate_aggr_lv) ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.alti <- felm(altg_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.slope <- felm(scale(slp_index) ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.access <- felm(scale(acc_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.workab <- felm(scale(wrk_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.nutri <- felm(scale(nut_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.water <- felm(water_area ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.veget <- felm(vegetation_area ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.cacao_i <- felm(cci_mean ~ p_coal| m_pair | 0 | m_pair, data=df_2000)
c00.cacao_r <- felm(ccr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.coffee_i <- felm(cfi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.coffee_r <- felm(cfr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.grass_i <- felm(pgi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.grass_r <- felm(pgr_mean ~ p_coal| m_pair | 0 | m_pair, data=df_2000)
c00.legum_i <- felm(pli_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.legum_r <- felm(plr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.maize_i <- felm(mzi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.maize_r <- felm(mzr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.rice_i <- felm(rci_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.rice_r <- felm(rcr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.soy_i <- felm(sbi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.soy_r <- felm(sbr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.sugar_i <- felm(sgi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.sugar_r <- felm(sgr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.amphib <- felm(amp_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.birds <- felm(brd_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.mammal <- felm(mml_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.urbdist <- felm(HubDist ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.urbexte <- felm(urb_extent ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.br230 <- felm(br230_distance_1993 ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.mainrds <- felm(main_rds_dist ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.roads <- felm(rds_dist ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.elf <- felm(elf ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.popc <- felm(pop95c_avg ~ p_coal | m_pair | 0 | m_pair, data=df_2000)
c00.popd <- felm(pop95d_avg ~ p_coal | m_pair | 0 | m_pair, data=df_2000)

# Balance - 2004
c04.ideol <- felm(ideo_imp ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.defo <- felm(deforested ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.rain <- felm(rain_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.eva <- felm(scale(eva_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.temp <-felm(temp_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.climaggr <- felm(scale(climate_aggr_lv) ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.alti <- felm(altg_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.slope <- felm(scale(slp_index) ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.access <- felm(scale(acc_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.workab <- felm(scale(wrk_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.nutri <- felm(scale(nut_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.water <- felm(water_area ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.veget <- felm(vegetation_area ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.cacao_i <- felm(cci_mean ~ p_coal| m_pair | 0 | m_pair, data=df_2004)
c04.cacao_r <- felm(ccr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.coffee_i <- felm(cfi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.coffee_r <- felm(cfr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.grass_i <- felm(pgi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.grass_r <- felm(pgr_mean ~ p_coal| m_pair | 0 | m_pair, data=df_2004)
c04.legum_i <- felm(pli_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.legum_r <- felm(plr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.maize_i <- felm(mzi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.maize_r <- felm(mzr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.rice_i <- felm(rci_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.rice_r <- felm(rcr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.soy_i <- felm(sbi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.soy_r <- felm(sbr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.sugar_i <- felm(sgi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.sugar_r <- felm(sgr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.amphib <- felm(amp_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.birds <- felm(brd_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.mammal <- felm(mml_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.urbdist <- felm(HubDist ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.urbexte <- felm(urb_extent ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.br230 <- felm(br230_distance_1993 ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.mainrds <- felm(main_rds_dist ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.roads <- felm(rds_dist ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.elf <- felm(elf ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.popc <- felm(pop95c_avg ~ p_coal | m_pair | 0 | m_pair, data=df_2004)
c04.popd <- felm(pop95d_avg ~ p_coal | m_pair | 0 | m_pair, data=df_2004)

# Balance - 2008
c08.ideol <- felm(ideo_imp ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.defo <- felm(deforested ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.rain <- felm(rain_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.eva <- felm(scale(eva_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.temp <-felm(temp_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.climaggr <- felm(scale(climate_aggr_lv) ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.alti <- felm(altg_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.slope <- felm(scale(slp_index) ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.access <- felm(scale(acc_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.workab <- felm(scale(wrk_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.nutri <- felm(scale(nut_mean) ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.water <- felm(water_area ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.veget <- felm(vegetation_area ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.cacao_i <- felm(cci_mean ~ p_coal| m_pair | 0 | m_pair, data=df_2008)
c08.cacao_r <- felm(ccr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.coffee_i <- felm(cfi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.coffee_r <- felm(cfr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.grass_i <- felm(pgi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.grass_r <- felm(pgr_mean ~ p_coal| m_pair | 0 | m_pair, data=df_2008)
c08.legum_i <- felm(pli_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.legum_r <- felm(plr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.maize_i <- felm(mzi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.maize_r <- felm(mzr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.rice_i <- felm(rci_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.rice_r <- felm(rcr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.soy_i <- felm(sbi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.soy_r <- felm(sbr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.sugar_i <- felm(sgi_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.sugar_r <- felm(sgr_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.amphib <- felm(amp_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.birds <- felm(brd_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.mammal <- felm(mml_mean ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.urbdist <- felm(HubDist ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.urbexte <- felm(urb_extent ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.br230 <- felm(br230_distance_1993 ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.mainrds <- felm(main_rds_dist ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.roads <- felm(rds_dist ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.elf <- felm(elf ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.popc <- felm(pop95c_avg ~ p_coal | m_pair | 0 | m_pair, data=df_2008)
c08.popd <- felm(pop95d_avg ~ p_coal | m_pair | 0 | m_pair, data=df_2008)


## Table A4: Moran’s I Test for Coalition Alignment on 40 Pre-Treatment Covariates in the 1996, 2000, 2004, and 2008 Mayoral Elections ----
# Define CRS and load grid cells shapefile
sad69 <- CRS("+proj=poly +lat_0=0 +lon_0=-54 +x_0=5000000 +y_0=10000000 +ellps=aust_SA +towgs84=-57,1,-41,0,0,0,0 +units=m +no_defs")
grid <- readOGR(dsn=path.expand("~/Dropbox/Research/Transamazonica/data/analysis/grid"), layer="hex_grid_new")
grid <- spTransform(grid, sad69)
# Merge with data frames by mayoral year
grid.96 <- subset(grid, ID %in% df_1996$ID) %>% merge(df_1996, by='ID', all.y=TRUE)
grid.00 <- subset(grid, ID %in% df_2000$ID) %>% merge(df_2000, by='ID', all.y=TRUE)
grid.04 <- subset(grid, ID %in% df_2004$ID) %>% merge(df_2004, by='ID', all.y=TRUE)
grid.08 <- subset(grid, ID %in% df_2008$ID) %>% merge(df_2008, by='ID', all.y=TRUE)
# Construct neighbors list from shapefile
reps <- 10
eps <- sqrt(.Machine$double.eps)
STRQ.96 <- system.time(for(i in 1:reps) a96 <- rgeos::gUnarySTRtreeQuery(grid.96)) / reps
STRQ.00 <- system.time(for(i in 1:reps) a00 <- rgeos::gUnarySTRtreeQuery(grid.00)) / reps
STRQ.04 <- system.time(for(i in 1:reps) a04 <- rgeos::gUnarySTRtreeQuery(grid.04)) / reps
STRQ.08 <- system.time(for(i in 1:reps) a08 <- rgeos::gUnarySTRtreeQuery(grid.08)) / reps
system.time(for(i in 1:reps) nb.96 <- poly2nb(grid.96, queen=TRUE, snap=eps, foundInBox=a96)) / reps + STRQ.96
system.time(for(i in 1:reps) nb.00 <- poly2nb(grid.00, queen=TRUE, snap=eps, foundInBox=a00)) / reps + STRQ.00
system.time(for(i in 1:reps) nb.04 <- poly2nb(grid.04, queen=TRUE, snap=eps, foundInBox=a04)) / reps + STRQ.04
system.time(for(i in 1:reps) nb.08 <- poly2nb(grid.08, queen=TRUE, snap=eps, foundInBox=a08)) / reps + STRQ.08
# Spatial weights for neighbors lists
wt.96 <- nb2listw(nb.96, zero.policy=TRUE) 
wt.00 <- nb2listw(nb.00, zero.policy=TRUE) 
wt.04 <- nb2listw(nb.04, zero.policy=TRUE) 
wt.08 <- nb2listw(nb.08, zero.policy=TRUE) 

# Moran I's - 1996
mi96.ideol <- moran.test(grid.96$ideo_imp, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi96.defo <- moran.test(grid.96$deforested, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.rain <- moran.test(grid.96$rain_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.eva <- moran.test(grid.96$eva_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.temp <- moran.test(grid.96$temp_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.rain <- moran.test(grid.96$temp_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.climaggr <- moran.test(grid.96$climate_aggr_lv, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.alti <- moran.test(grid.96$altg_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.slope <- moran.test(grid.96$slp_index, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.access <- moran.test(grid.96$acc_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.workab <- moran.test(grid.96$wrk_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.nutri <- moran.test(grid.96$nut_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.water <- moran.test(grid.96$water_area, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.veget <- moran.test(grid.96$vegetation_area, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi96.cacao_i <- moran.test(grid.96$cci_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.cacao_r <- moran.test(grid.96$ccr_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.coffee_i <- moran.test(grid.96$cfi_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.coffee_r <- moran.test(grid.96$cfr_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.grass_i <- moran.test(grid.96$pgi_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.grass_r <- moran.test(grid.96$pgr_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.legum_i <- moran.test(grid.96$pli_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.legum_r <- moran.test(grid.96$plr_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.maize_i <- moran.test(grid.96$mzi_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.maize_r <- moran.test(grid.96$mzr_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.rice_i <- moran.test(grid.96$rci_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.rice_r <- moran.test(grid.96$rcr_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.soy_i <- moran.test(grid.96$sbi_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.soy_r <- moran.test(grid.96$sbr_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.sugar_i <- moran.test(grid.96$sgi_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.sugar_r <- moran.test(grid.96$sgr_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.amphib <- moran.test(grid.96$amp_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.birds <- moran.test(grid.96$brd_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.mammal <- moran.test(grid.96$mml_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.amphib <- moran.test(grid.96$mml_mean, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.urbdist <- moran.test(grid.96$HubDist, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.urbexte <- moran.test(grid.96$urb_extent, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi96.br230 <- moran.test(grid.96$br230_distance_1993, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.mainrds <- moran.test(grid.96$main_rds_dist, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.roads <- moran.test(grid.96$rds_dist, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.elf <- moran.test(grid.96$elf, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi96.popc <- moran.test(grid.96$pop95c_avg, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi96.popd <- moran.test(grid.96$pop95d_avg, wt.96, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")

# Moran I's - 2000
mi00.ideol <- moran.test(grid.00$ideo_imp, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi00.defo <- moran.test(grid.00$deforested, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.rain <- moran.test(grid.00$rain_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.eva <- moran.test(grid.00$eva_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.temp <- moran.test(grid.00$temp_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.rain <- moran.test(grid.00$temp_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.climaggr <- moran.test(grid.00$climate_aggr_lv, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.alti <- moran.test(grid.00$altg_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.slope <- moran.test(grid.00$slp_index, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.access <- moran.test(grid.00$acc_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.workab <- moran.test(grid.00$wrk_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.nutri <- moran.test(grid.00$nut_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.water <- moran.test(grid.00$water_area, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.veget <- moran.test(grid.00$vegetation_area, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi00.cacao_i <- moran.test(grid.00$cci_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.cacao_r <- moran.test(grid.00$ccr_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.coffee_i <- moran.test(grid.00$cfi_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.coffee_r <- moran.test(grid.00$cfr_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.grass_i <- moran.test(grid.00$pgi_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.grass_r <- moran.test(grid.00$pgr_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.legum_i <- moran.test(grid.00$pli_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.legum_r <- moran.test(grid.00$plr_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.maize_i <- moran.test(grid.00$mzi_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.maize_r <- moran.test(grid.00$mzr_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.rice_i <- moran.test(grid.00$rci_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.rice_r <- moran.test(grid.00$rcr_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.soy_i <- moran.test(grid.00$sbi_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.soy_r <- moran.test(grid.00$sbr_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.sugar_i <- moran.test(grid.00$sgi_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.sugar_r <- moran.test(grid.00$sgr_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.amphib <- moran.test(grid.00$amp_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.birds <- moran.test(grid.00$brd_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.mammal <- moran.test(grid.00$mml_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.amphib <- moran.test(grid.00$mml_mean, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.urbdist <- moran.test(grid.00$HubDist, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.urbexte <- moran.test(grid.00$urb_extent, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi00.br230 <- moran.test(grid.00$br230_distance_1993, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.mainrds <- moran.test(grid.00$main_rds_dist, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.roads <- moran.test(grid.00$rds_dist, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.elf <- moran.test(grid.00$elf, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi00.popc <- moran.test(grid.00$pop95c_avg, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi00.popd <- moran.test(grid.00$pop95d_avg, wt.00, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")

# Moran I's - 2004
mi04.ideol <- moran.test(grid.04$ideo_imp, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi04.defo <- moran.test(grid.04$deforested, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.rain <- moran.test(grid.04$rain_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.eva <- moran.test(grid.04$eva_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.temp <- moran.test(grid.04$temp_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.rain <- moran.test(grid.04$temp_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.climaggr <- moran.test(grid.04$climate_aggr_lv, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.alti <- moran.test(grid.04$altg_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.slope <- moran.test(grid.04$slp_index, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.access <- moran.test(grid.04$acc_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.workab <- moran.test(grid.04$wrk_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.nutri <- moran.test(grid.04$nut_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.water <- moran.test(grid.04$water_area, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.veget <- moran.test(grid.04$vegetation_area, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi04.cacao_i <- moran.test(grid.04$cci_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.cacao_r <- moran.test(grid.04$ccr_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.coffee_i <- moran.test(grid.04$cfi_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.coffee_r <- moran.test(grid.04$cfr_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.grass_i <- moran.test(grid.04$pgi_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.grass_r <- moran.test(grid.04$pgr_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.legum_i <- moran.test(grid.04$pli_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.legum_r <- moran.test(grid.04$plr_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.maize_i <- moran.test(grid.04$mzi_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.maize_r <- moran.test(grid.04$mzr_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.rice_i <- moran.test(grid.04$rci_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.rice_r <- moran.test(grid.04$rcr_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.soy_i <- moran.test(grid.04$sbi_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.soy_r <- moran.test(grid.04$sbr_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.sugar_i <- moran.test(grid.04$sgi_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.sugar_r <- moran.test(grid.04$sgr_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.amphib <- moran.test(grid.04$amp_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.birds <- moran.test(grid.04$brd_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.mammal <- moran.test(grid.04$mml_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.amphib <- moran.test(grid.04$mml_mean, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.urbdist <- moran.test(grid.04$HubDist, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.urbexte <- moran.test(grid.04$urb_extent, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi04.br230 <- moran.test(grid.04$br230_distance_1993, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.mainrds <- moran.test(grid.04$main_rds_dist, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.roads <- moran.test(grid.04$rds_dist, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.elf <- moran.test(grid.04$elf, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi04.popc <- moran.test(grid.04$pop95c_avg, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi04.popd <- moran.test(grid.04$pop95d_avg, wt.04, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")

# Moran I's - 2008
mi08.ideol <- moran.test(grid.08$ideo_imp, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi08.defo <- moran.test(grid.08$deforested, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.rain <- moran.test(grid.08$rain_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.eva <- moran.test(grid.08$eva_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.temp <- moran.test(grid.08$temp_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.rain <- moran.test(grid.08$temp_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.climaggr <- moran.test(grid.08$climate_aggr_lv, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.alti <- moran.test(grid.08$altg_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.slope <- moran.test(grid.08$slp_index, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.access <- moran.test(grid.08$acc_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.workab <- moran.test(grid.08$wrk_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.nutri <- moran.test(grid.08$nut_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.water <- moran.test(grid.08$water_area, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.veget <- moran.test(grid.08$vegetation_area, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi08.cacao_i <- moran.test(grid.08$cci_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.cacao_r <- moran.test(grid.08$ccr_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.coffee_i <- moran.test(grid.08$cfi_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.coffee_r <- moran.test(grid.08$cfr_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.grass_i <- moran.test(grid.08$pgi_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.grass_r <- moran.test(grid.08$pgr_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.legum_i <- moran.test(grid.08$pli_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.legum_r <- moran.test(grid.08$plr_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.maize_i <- moran.test(grid.08$mzi_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.maize_r <- moran.test(grid.08$mzr_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.rice_i <- moran.test(grid.08$rci_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.rice_r <- moran.test(grid.08$rcr_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.soy_i <- moran.test(grid.08$sbi_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.soy_r <- moran.test(grid.08$sbr_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.sugar_i <- moran.test(grid.08$sgi_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.sugar_r <- moran.test(grid.08$sgr_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.amphib <- moran.test(grid.08$amp_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.birds <- moran.test(grid.08$brd_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.mammal <- moran.test(grid.08$mml_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.amphib <- moran.test(grid.08$mml_mean, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.urbdist <- moran.test(grid.08$HubDist, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.urbexte <- moran.test(grid.08$urb_extent, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi08.br230 <- moran.test(grid.08$br230_distance_1993, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.mainrds <- moran.test(grid.08$main_rds_dist, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.roads <- moran.test(grid.08$rds_dist, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.elf <- moran.test(grid.08$elf, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided", na.action=na.omit)
mi08.popc <- moran.test(grid.08$pop95c_avg, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")
mi08.popd <- moran.test(grid.08$pop95d_avg, wt.08, zero.policy=TRUE, randomisation=TRUE, alternative="two.sided")

